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SUMMARY 

The fitting of mathematical models to physiological systems can be tedious and difficult, whether 
one uses analog or digital computer methods. Both methods have their pros and cons depending 
on the available hardware and software and on the type of modeling. In recent years many digital 
simulation languages have been written combining analog-like and digital features to facilitate 
modeling, but, for a variety of reasons, none of these was suitable for our applications. We, 
therefore, designed a new digital simulation control system, SIMCON, which is described in this 
paper. 

The primary objectives were to provide: 

1. Maximum man-machine interaction at run-time, including visual displays, digital control, 
and both continuous analog and digital parameter adjustment 

2. The ability to generate solutions and to fit them to experimental data or other theoretical 
curves with a minimum of computer memory 

3. The option to use a mathematically oriented language, FORTRAN, and block operators 
with variable input/output. 

The result is a relatively general and simple simulation system which is easy to use and has wide 
versatility. 


INTRODUCTION 

The fitting of mathematical models to physiological systems can be tedious or difficult when 
using either analog or digital methods. Typically, analog simulations require relatively long 
setup time and some care in time and amplitude scaling, but they provide high solution 
speed, visual output, and easy changing of parameters. While digital solutions are easy to 
program, and while hard copy output and greater accuracy may be obtained digitally, the 
serial nature of digital computation makes solutions slow, and there are seldom adequate 
means for making immediate visual assessments of the solutions or for making adjustments 
to the model or its parameters. Many digital simulation languages that combine features of 
both analog and digital machines have been written in recent years in attempts to facilitate 
the solution of such problems as model fitting. However, after studying several of these 
languages, we concluded that none was suitable for our applications, for a variety of reasons 
—most lacked on-line operational flexibility, many were not compatible with our particular 
hardware and software configurations, some had extravagant memory requirements, 
involved long or special setup procedures, lacked a variety of input/output characteristics, 
etc. 

The objectives of the digital simulation system (SIMCON) were (1) maximum man-machine 
interaction at run time including visual displays, digital control, and both continuous analog 
and digital parameter adjustment, (2) the ability to generate solutions and to fit them to 
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experimental data or other theoretical curves without using large amounts of computer 
memory, and (3) the combined use of a mathematically oriented language, FORTRAN, and 
block operators with variable inputs and outputs. 

DESCRIPTION 

General 

SIMCON was designed for use with a CDC 3200 computer, a time-sharing monitor (a 
modified MEDLAB1), one tape drive, one disk pack, a printer, a card reader, analog-to- 
digital and digital-to-analog converters, and a Calcomp plotter. Program execution is 
initiated from a peripheral computer station consisting of a decimal keyboard and a memory 
oscilloscope (Figure 1). SIMCON, written in basic assembly language, BAP, provides a 
means of investigator-model interaction via digital and analog input/output, while taking full 
advantage of the time-sharing monitor. 

The mathematical model is programmed in FORTRAN II and after compilation is stored on 
the magnetic disk with an identifying number. When the FORTRAN program number is 
entered with the decimal keyboard, SIMCON loads the program into memory and controls 
its execution. Control parameters and data that must be addressed by both SIMCON and the 
FORTRAN program are stored in the common area of memory. 

The control of SIMCON (see flow chart. Figure 2) is almost entirely accomplished by use of 
parameter options, i.e., the investigator manipulates the flow of the program by setting 
certain controlling parameters to the “on” or “off’ position. Some of these on/off options are 
illustrated as decision blocks (black diamonds) accompanied by the controlling parameter 
number. One hundred “parameters” (Pioo to P 200 ) are associated with control functions 
within SIMCON and the first hundred are available for use in the user’s FORTRAN 
program (see Figure 3). Any of the parameters can be displayed on the memory oscilloscope 
or written on magnetic tape or the Calcomp plotter. 

The control parameters and the user’s parameters can be adjusted either prior to or during a 
solution, by any of several means: 

1. Card input prior to a solution 

2. Decimal keyboard in conjunction with a memory oscilloscope at one of the remote 
stations 

3. Analog potentiometers 

4. Automatic parameter adjustment loops 

5. The FORTRAN program itself. 

Figures 4a and 4b are typical messages on the memory oscilloscope showing the means of 
computer-to-user communication afforded by SIMCON. The combination of digits shown 
on the left of Figure 4a is entered on the keyboard at the remote terminal to direct the flow 
of the program. An entry on the keyboard preceded by a 1 is recognized as a parameter 
index; an entry preceded by a 2 is recognized as a new parameter value. If the user enters an 
incorrect value from the keyboard, he may depress the SOM (start of message) key, which 
allows the entry of a new value. Figure 4b shows the oscilloscope message received when 
parameter 57 was set to a new value from the remote terminal. 

Modeling program 

The program for mathematical modeling is written in FORTRAN with all the standard 
features (I/O, function subroutines, etc.) available, but control is exerted by the user at a 
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peripheral station. Subroutines may be written in the usual way, and a growing number of 
block operator subroutines (integrators of several types, linear and parabolic interpolation, 
fixed and continuously variable delays, function generators, etc.) are available as 
FORTRAN subroutines. While we have chosen to use FORTRAN as the working language, 
there is no reason that we could not have used another, more specialized simulation 
language. 

The FORTRAN program is similar to those used with the SCi continuous system simulation 
language 2 in that it has three parts: (1) an initial region, (2) a dynamic or solution-generating 
region, and (3) a terminal region. 

The initial region is executed once at the beginning of a solution and sets the initial 
conditions for the solution. If variable, these initial conditions can be entered with the 
decimal keyboard. 

The terminal region of the FORTRAN program is used to make a quantitative assessment of 
the solution and to produce a permanent record of the solution on chosen output devices. 

Any variable that is selected at run time can also be displayed graphically on the memory 
oscilloscope during the solution, with new points being added with each step increase of the 
independent variable (usually time). Parameters used in the FORTRAN program can be 
altered via the keyboard prior to and during generation of the solution to determine solution 
sensitivities to parameter changes. Experimental or other data arrays that have been 
previously stored on the magnetic disk may be displayed on the memory oscilloscope for 
visual comparison with the solutions obtained from a model. The stored arrays may also be 
variable functions used in obtaining the solution, usually after interpolation. When repetitive 
solutions are being obtained, up to 10 computed values for each solution may be written on 
digital tape for later display or analysis. (The number 10 was arbitrarily chosen to meet our 
present requirements and could readily be increased.) 

Control loops 

External to and nested around the innermost solution loop are three other loops (see Figure 
2), each of which may be functional or nonfunctional as determined by parameter options. 
The first is an “analog potentiometer loop” by means of which parameters take values either 
prior to or during a solution from 4 potentiometers whose outputs are converted to digital 
values for use during a solution. The analog signals (0 to 3 volts DC) are generated by 
batteries within the box shown above the oscilloscope in Figure 1 and are transmitted to the 
digital computer via a 10-bit analog-to-digital converter. Each potentiometer has three 
controlling parameters associated with it. The first specifies which one of the two hundred 
parameters is to be changed during a solution, and the other two are chosen maximum and 
minimum values for the parameter which is to be altered. The value of the specified 
parameter is the chosen minimum plus the fraction of the potentiometer range times the 
difference between the maximum and the minimum. The potentiometers allow a more rapid 
adjustment of parameter values between iterations than could be obtained by digital entry of 
new values, and they can be continuously varied during a solution. 

The other two outer loops are used to make predetermined arithmetic or logarithmic 
adjustments of 2 selected parameters. Parameters may also be adjusted under program 
control—for example, after assessment of the solution in the third portion of the user’s 
program. 
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Output 

Hard copy graphical output for display of variables is obtainable by means of (a) the 
memory scope, which can be photographed, (b) Calcomp graphs, or (c) printer plots. The 
same methods are available for plotting results from many solutions stored on digital tape— 
for example, for plotting the amplitude of output of an operator against the frequency of a 
repetitive input function. The arrays stored on tape may also be displayed rapidly on the 
memory scope using arithmetic or logarithmic scaling, thus allowing a quick check of 
results or choice of appropriate scaling for permanent graphing. 


Applications 

Many processes in biological research can be approached by describing the transfer function 
of the system whose input and output can be measured experimentally. This is the classical 
problem of system identification; it may be attacked in various ways, depending in part on 
whether the system is a simple open-loop network or a multiple-feedback closed-loop 
system. An illustration is the problem of obtaining a description of the transport function, 
the probability density function of transit times of material carried in the bloodstream 
between an upstream point A and a downstream point B. Let some indicator which travels 
with the blood be injected upstream and its concentration in the blood at A and B be 
recorded as Aft) and B(t), respectively. Using SIMCON, the transport function from A to B 
can be obtained accurately and with relative ease based on the assumption that this transport 
function can be approximated with what may be described as a variable-order low-pass filter 
whose impulse response is overdamped. This last requirement is important since negative 
concentrations cannot exist. 

Figure 5 is a block diagram of the FORTRAN program for the variable-order model used to 
simulate the transport function. Aft) and Bft) represent the digitized upstream and 
downstream indicator dilution curves that have been read into memory from the magnetic 
disk. Many pairs of curves may reside on the disk, and a particular pair is selected with 
SIMCON parameter options. Operators FIRST and SECORD approximate dispersion plus 
transport delay, while DELAY is an operator that approximates delay alone. All operator 
parameters, initial conditions, damping coefficients (ZETA), natural frequencies (OMEGA), 
time constants (TAU), and delays (DLAY) are set according to SIMCON parameter values. 
Parameters 89 and 90 are the outputs of scaling summers, and parameter 71 is the “integral” 
of the square of the differences between two functions used to measure curve ’’fit.” The 
linkages between these operators are quite variable, and as long as the step size in the 
independent variable (time) is small, then the order of listing in the program is of little 
consequence. A linkage such as given by the dashed lines must be used, certain operators 
remaining unused. When the parameters of the model have been adjusted so that its output 
response to an input having the form of the upstream dye curve closely approximates the 
recorded downstream dye curve, then the impulse response of the model can be generated by 
substituting a short duration pulse, parameter 82, as the input to the model. 

The FORTRAN model program is shown below to illustrate its simplicity. The initializing, 
solution generating, and solution assessing sections of the program are shown, as is the 
scheme for linking the operators. 
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PROGRAM 75«0 

COMMON I A, 1B 

DIMENSION I A(7 <5 ) , 1B(745),P(200) 

p is the array df variable control parameters. 

T IME*Pi130 ) 

TIME IS THE INDEPENDENT VARIABLE OF THIS SOLUTION, 

IF(TIHE) 1.1*2 

....PARAMETER INITIALIZATION SECTION AT TIME. 0.0 ONLV.. 

TIME CONSTANTS FOR LAO OPERATORS, 


1 ZE T Al.P(11) 
ZETA2«P(2l> 
ZETA3«P(3l) 


OMEGA1.PI12) 
OMEGA2.PI22) 
0HEGA3.P132 > 
TAUl.PIAl) 
TAU2*PM6) 
TAU3.PU1) 

DL AY1«P{56 > 
DLAY2.PI61) 
DLAr3oP(66) 

YZEROl.P113) 
YZER02'PI23) 
YZER03*P<33> 
YZER04.P(A3) 
VZER05«P(4eI 
YZER06.P <33) 
DYZER01.P<14> 
DYZER02.P(24} 
OY ZERD3.PI34 > 


NS1• P(91> 
NS2« P192 > 
NS3« P(93 > 
NS4. P(94> 

Nl.PIlI 

N2.P<2) 

N3*P<3l 

N4.PC4) 

n5«P<5) 

N6.PI6) 

N7 *P(71 
N8.PI8I 


(SEC) FOR DELAY OPERATORS. 


INITIAL CONDITIONS FOR DIFFERENTIAL OPERATORS. 


N9. 


(9) 


PULSE • P(15) 
PULOUR.PIie) 


INDICES OF VARIABLES USED AS INPUTS TO OPERATORS. 


INPUT PULSE AMPLITUDE 
INPUT PULSE DURATION, 

INDICES OF VARIABLES TO BE COmPaReo. 


NFNCT.P(17) 

KOMPARE « P(18) 

SUMSO.O,0 

........SOLUTI OK GENERATION SECT ION«.«••••••.••••••••••••••••••••••••••..... 

INTERP INTERPOLATES DaTa ARRAYS FROM 01SK 
(IA AND 16) TO FINO VALUES *T PROPER INTERVALS, 

2 INPUT.INTERPCIA) 

OUTPUT*INTERPI|B> 

P ( 80 > * I NP.UT 
P(81)sOUTPLT 

P(82).PULSE 

IF(TIME-PULDLR) 22,22,21 

21 P (82 > *0. 

SECORD IS a SECOND ORDER OPERaTdR, 

22 P<10).SECORD(P(N1).ZETA1,0MEGA1.YZER01,DYZERD1.1) 
P(20).SECORO(P<N2),ZETa2,OMEGa2,TZERO2 i DYZER02.2) 
P<30).SECORD<P(N3).ZETa3.DHEGa3,YZERO3,DYZERO3,3) 

FIRST IS A LAG OPERATOR. 

P(40) .FIRST(P(N4),TAU1.YZER04,1) 

P(43).F|RST(P<N5),TAU2.YZER05,2) 

PI 50)=FIRSI(P(N6),Tau3.YZER06,3> 

DELAY |S A PURE, FIVEO DELAY OPERATOR, 

PI 55 ) * DEL*Y«P(N7), 0LAY1.PI55), 1) 

P(60)« DELAY(P(N6), DLAY2, P(60),2) 

PI65)*DELA»(P(N9), DLAY3, P C 65 >,3 > 

P(69) ANO P < 90) ARE OUTPUTS DF SUMMERS. 

P < 89 ) * PI 95)*P(NSl1 . P(96)«P(NS2) 

P(90 ) . PI 97)•PINS3> * PI98).P(NS4) 

SUMSO IS The TEST STATISTIC OP THE FIT. 

SUMSO. SUMSO . IP(KDMPARE)-PINFNCT)).*2 
P< 71).SUMSC 

P(129) IS THE TOTAL SOLUTION TIME, 

IFIPI129)-TIME) 3,3,4 

........SOLUTION ASSESSMENT SECT ION--LAST PASS ONLY........«••' 

3 C0EF«SUmSQ/<T1mE/PI131)) 

PI131) IS THE TIME INCREMENT, 

P(7 2)sCOEF 

IF P<19 9) IS POSITIVE. PRINT P ARRAY, 

IF <P<199)) 4.4,5 
5 NP.PI199) 

PRINT 199. IP(| ) ,1*1,NP) 

199 FORMATC1X.E9.3.9E10.3) 

4 STOP 
ENO 


Figure 6 shows an upstream indicator dilution curve, and an output of the model, 

labeled C^t) * hit), which was compared with an experimentally recorded downstream 
indicator dilution curve The parameters of the model were varied, using the keyboard 
at a remote computer terminal, until the differences between the theoretical and the 
experimental output curves were minimized; when the differences were less than a 
predetermined tolerance, the model was said to approximate the transport function, h(t). A 
third-order filter plus a time delay was used to obtain this fit, one SECORD (ZETA = 0.85, 
OMEGA = 1.20), one FIRST (TAU = 0.70), and one DELAY (DLAY = 5.20) all in series. 

The goodness of fit between C#(7) and C^t) * hit) was determined visually on the memory 
oscilloscope and improved by minimizing the sums of squares of the differences between 
the computed solution and the experimental curve. This particular transport function 
description required approximately five minutes for parameter adjustment. Since the updated 
parameter list is stored with the FORTRAN program on the magnetic disk, successive fits of 
similar functions require less iteration time for the selection of an optimum model and 
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model parameters. This program thus contained more operators than were necessary, but this 
excess provided wide latitude in choosing and testing a variety of models. Having chosen a 
suitable model as a result of such testing, one should then rewrite the program omitting the 
unnecessary operators. 

In feedback systems such as might be formulated at least crudely using a general program 
(such as program 7540 or Figure 5), the phase lags introduced by inappropriate sequencing 
of the operators can be partly overcome by using a very small step size in obtaining the 
solutions. After choosing a better defined model the sequencing of the operators should be 
changed to minimize undesired phase lags and to permit the use of larger step sizes. 

In another application, SIMCON was used to define magnitudes of potential errors in 
estimating flow, mean transit time, dispersion, and volume by the standard dye dilution 
techniques introduced by fluctuations in flow. 3 The adjustable parameters included the 
frequency, phase and amplitude of the sinusoidal flow, shaping parameters of the dispersive 
delay operator, the average flow, the form of the simulated tracer input, and many others 
including the selection of output variables and output devices. Other applications have been 
in testing models for capillary-tissue exchange and analyzing data obtained from dogs for 
permeability of capillary membranes in the heart 4 and for the extent of counter-current 
diffusional exchange of water in the heart. 5 

CONCLUSION 

In general, SIMCON is very useful in finding the solution to any problem where repeated 
man-machine interaction is essential because logical adjustments are difficult or impossible 
to define in advance. Most obvious are the iterative methods of solution where error 
functions are difficult to program or define before the solution is initiated. SIMCON 
provides the programming ease, accuracy, memory, and hard copy output of digital 
computing and yet allows the man-machine interaction and run-time flexibility usually 
associated with analog computing. 
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Figure 1. 

A computer terminal consisting of an oscilloscope, a keyboard and a potentiometer- 
controlled parameter entry device 
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Figure 2. 

Block diagram of SIMCON. The numbers beside the decision blocks (black diamonds) 
indicate which parameter is set to control the program route. For example, P 109 would be set 
positive if the user wishes to read the analog line voltages controlled by the analog 
potentiometers. 
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Figure 3. 

SIMCON parameter functions 
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EXIT 

PARAM INDEX I 
ENTER PARAM I 


Figure 4. 

(a) One of the messages transmitted by SIMCON prior to an operator decision. SOM and 
EOM are keyboard codes for start and end of message, respectively, (b) Reply from 
SIMCON after a new value has been entered for parameter 57. 
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BLOCK DIAGRAM OF PROGRAM 7540 



Figure 5. 

Diagram of operators available in program 7540. Dashed lines show possible configuration 
chosen at run time. 
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CIRCULATORY TRANSPORT FUNCTION 
approximated with a third order model plus delay 



TIME (seconds) 


Figure 6. 

Upstream, C 1 \{t), and downstream, Cg(f), dilution curves are shown along with a theoretical 
downstream curve (dashed line) obtained using a simulated transport function, h(t) 
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